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ABSTRACT 

We present a new method for determining the local dark matter density using kinematic data 
for a population of tracer stars. The Jeans equation in the ^-direction is integrated to yield 
an equation that gives the velocity dispersion as a function of the total mass density, tracer 
density, and the ‘tilt’ term that describes the coupling of vertical and radial motions. We 
then fit a dark matter mass profile to tracer density and velocity dispersion data to derive 
credible regions on the vertical dark matter density profile. Our method avoids numerical 
differentiation, leading to lower numerical noise, and is able to deal with the tilt term while 
remaining one dimensional. In this study we present the method and perform initial tests 
on idealised mock data. We also demonstrate the importance of dealing with the tilt term 
for tracers that sample > 1 kpc above the disc plane. If ignored, this results in a systematic 
underestimation of the dark matter density. 
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1 INTRODUCTION 

Dark matter (DM) is an elusive component of the cosmos. Its exis¬ 
tence has long been inferred from its gravitational interactions with 
ordinary matter on scales ranging from dwarf galaxies and galactic 
rotation curves to lensing galaxies, galaxy clusters, and the Uni¬ 
verse as a whole (for reviews see e.g. Rubin 1993; Jungman et al. 
1996; Bergstrom 2000; Bertone et al. 2005; Bertone 2010; Massey 
et al. 2010). However its exact nature remains unknown. 

Understanding the distribution of DM in cosmological struc¬ 
tures, and in particular in the Milky Way, is of great importance to 
testing the standard A-Cold Dark Matter (A-CDM) cosmological 
model; to make predictions - and allow the interpretation - of ex¬ 
periments that seek to detect DM; and to probe more exotic models 
for DM (e.g. Read 2014). Here we will focus our attention on the 
local DM density - i.e. the average density of DM over a small vol¬ 
ume (typically of size ~ lOOpc - lOOOpc) centred around the Sun. 
The local DM density is a vital input for so-called direct DM exper¬ 
iments, based on the detection of the recoil energy of nuclei struck 
by DM particles streaming through underground detectors, as the 
expected event rate is proportional to the local DM density and the 
DM-nucleon cross section. In addition, the local DM density is im¬ 
portant for indirect searches that look for neutrinos produced by 
the annihilation of DM particles trapped at the centre of Sun (see 
e.g. Klasen et al. 2015 and references therein for a recent update on 
direct and indirect searches). The results of these searches are also 
used in explorations of theoretical parameter space, such as those of 


Supersymmetry (e.g. Allanach & Lester 2006; Buchmueller et al. 
2014; Strege et al. 2014; Roszkowski et al. 2015; Bertone et al. 
2016). 

Methods of determining the local DM density can be divided 
into two categories: those utilising measurements of stars in a small 
volume around the Sun (e.g. Kapteyn 1922; Jeans 1922; Oort 1932; 
Bahcall 1984; Kuijken & Gilmore 1989b, 1991; Creze et al. 1998; 
Garbari et al. 2012; Bovy & Tremaine 2012; Smith et al. 2012; 
Zhang et al. 2013; Bienayme et al. 2014) and those utilising a host 
of dynamical tracers - particularly the rotation curve - to constrain 
more global mass models of the Milky Way (Dehnen & Binney 
1998; Weber & de Boer 2010; Catena & Ullio 2010; Salucci et al. 
2010; McMillan 201 1; Nesti & Salucci 2013; Piffl et al. 2014; Iocco 
et al. 2015; Pato & Iocco 2015; Pato et al. 2015). Following the no¬ 
tation of Read (2014) we denote the local DM density derived from 
local measurements by Pdm> an( I those extrapolated from rotation 
curves assuming a spherical DM halo as poM.ext- With the advent of 
large survey data, local and global methods are beginning to con¬ 
verge (Bovy & Rix 2013; Piffl et al. 2014). 

The comparison of p DM to p DMext can provide insight into the 
shape of the Milky Way’s DM halo, and thus into the formation of 
the galaxy (Read 2014). If p DM from local measurements is smaller 
than that extrapolated from global measurements, poM.ext. then this 
implies a prolate halo, while the opposite would imply a squashed, 
oblate halo. While the former is produced by DM only N-body sim¬ 
ulations, the introduction of baryons into the simulations produces 
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the latter (Katz & Gunn 1991; Dubinski 1994; Debattista et al. 
2008; Read et al. 2009). 

If p DM is derived from tracers distributed vertically in the plane 
it is possible to probe not only the local DM density but also its ver¬ 
tical profile, PdmOi)- This could potentially reveal deviations from 
a spherical halo profile, which can form in a number of ways. First, 
the DM halo should respond to the formation of the baryonic disc 
by flattening into an oblate halo (Katz & Gunn 1991; Dubinski 
1994; Debattista et al. 2008; Read et al. 2009). Second, the accre¬ 
tion of subhalos after the formation of the baryonic disc can also 
lead to the formation of a ‘dark disc’ that co-rotates with the bary¬ 
onic disc (Lake 1989; Read et al. 2008, 2009; Purcell et al. 2009; 
Ling et al. 2010; Pillepich et al. 2014). Gravitationally these two 
features, an oblate halo and a dark disc, are indistinguishable and 
degenerate. One method to distinguish the two is to hunt for the 
chemically and kinematically distinct sub-halo stars that would ac¬ 
company the accreted dark disc, but not the contracted DM halo. 
Recent Gaia-ESO observational data find no evidence for such stars 
(Ruchti et al. 2014; Ruchti et al. 2015), suggesting that the Milky 
Way has had a rather quiet past since its disc formed, and therefore 
does not host a significant accreted dark disc. If correct, then any 
gravitationally detected non-sphericity must imply a non-rotating 
locally oblate halo. An enhancement in pdm over poM.ext is also a 
feature in some more exotic models, such as Partially Interacting 
Dark Matter (PIDM) (Fan et al. 2013), and some modified theories 
of gravity (e.g. Read & Moore 2005; Nipoti et al. 2007). 

Measurement of the local DM density dates back almost 100 
years (Kapteyn 1922). The history of this measurement is one of 
both increasingly precise data and decreasingly strong assump¬ 
tions. As such, the error bars on p DM have not always shrunk with 
time, as better data often allows one to discard previous assump¬ 
tions. With the advent of Gaia data from 2016 onwards we will 
have access to high precision data on individual stars, and the pri¬ 
mary uncertainty in the determination of pdm will be systematic 
model uncertainties. Recent results for the local DM density in¬ 
clude Bienayme et al. (2014), who used RAVE data to derive a 
value of pdm = 0.0143 ±0.001 lM 0 pc~ 3 = 0.542 ± 0.042 GeVcnT 3 , 
and Zhang et al. (2013) who used SDSS/SEGUE data to find 
Pdm = 0.0065 ± 0.0023M o pc -3 = 0.25 ± 0.09GeVcm~ 3 . Note that 
these two results do not overlap within their stated uncertainties. 
The significance of this discrepancy is hard to interpret though, 
as they each use different data sets and analysis techniques, and 
both methodologies rely on rather different assumptions. To make 
progress we should endeavour to minimise the number of assump¬ 
tions made, and apply the same analysis techniques to both data 
sets. 

In this paper we make progress towards the first goal of re¬ 
ducing model assumptions. To achieve this we introduce a one¬ 
dimensional Jeans analysis method to probe the local DM density 
using the vertical motions of tracer stars. We construct a representa¬ 
tion of the tracer density v, and also allow for a DM density profile 
that is more complex than simply constant with height as previous 
local studies have assumed. Additionally, we deal with the so-called 
tilt term, which links radial and vertical motions of the tracer stars. 
This term is crucial to understand stellar motions at high-z, where 
the baryonic contribution falls off and DM becomes increasingly 
important. Using the vertical Jeans equation we calculate the ve¬ 
locity dispersion cr.(z) for each mass model, and then fit to data in 
v, cr., and cr Rz using MultiNest (Feroz & Hobson 2008; Feroz et al. 
2009; Feroz et al. 2013). We test this method on mock data sets, and 
explore the impact of tracer star sample size, the tilt term, and non¬ 
constant DM profiles. The mock data for this paper are ‘as good as 


it gets’ - we assume the population to have no observational biases 
and there to be no measurement error on individual stars - allowing 
us to isolate the effects of sampling error and model uncertainties. 
This is similar to what we will have with Gaia data, where the mea¬ 
surement errors on stars will be small compared to sampling error. 
We will explore the effect of adding realistic Gaia uncertainties to 
our method in a separate work. 

We reiterate that this method is one dimensional, allowing us 
to keep assumptions to a minimum. However, we show that we are 
still able to deal with tilt and high-z data, which usually requires 
a 2-dimensional method. The key disadvantages of our method are 
that we bin data and thus lose information, and that our strictly local 
method cannot take advantage of the Tong lever arm’ of a global 
model that would ensure, for example, that the DM density in ra¬ 
dial slices close to the ‘local volume’ is continuous and smooth. 
However the effect of these disadvantages is to overestimate the er¬ 
rors on pdm. which is acceptable as we aim for a conservative and 
robust estimate of pdm- 

The paper is organised as follows: in Section 2 we introduce 
our method, covering the Jeans equation based mathematical for¬ 
malism, our treatment of the rotation curve and tilt term, our our 
descriptions of the elements of our mass and tracer density models. 
We then outline our statistical analysis, introducing the framework 
for Bayesian parameter estimation and the MultiNest nested sam¬ 
pling code. In Section 3 we describe the array of mock data sets 
we use to test out methods. In Section 4 we present and discuss 
the results of our analyses, before finally concluding the paper with 
Section 5. 


2 METHOD 

The broad picture of our problem is this: we have quantities derived 
from the motions of the tracer stars, namely v, the tracer density, 
cr z , the vertical velocity dispersion, and cr Rz , the (R, z) velocity dis¬ 
persion. Then we have elements of the mass profile, one of which 
is unknown - the DM density Pdm. and the other which is known 
within a band of uncertainty - the baryon density p baryon . In this sec¬ 
tion we first derive the equations to link these quantities, and then 
describe how each is modelled. 


2.1 Deriving a general ID Jeans method 

A population of ‘tracer’ stars moving vertically near the Sun obeys 
the collisionless Boltzmann equation: 

df df 

dt = ^ + Vx/-v-V v /-V,® = 0 (1) 

where /(x, v) is the stellar distribution function; x and v are the 
positions and velocities, respectively; and ® is the gravitational po¬ 
tential. 

If the system is in dynamic equilibrium (steady state), then we 
may neglect the partial time derivative of /. Thus for equilibrium 
tracers, if we know their phase space distribution function /, then 
we can directly measure the gravitational force from equation (1). 
In practice, however, this is hard because / is six-dimensional (even 
a million stars gives only 10 sample points per dimension) and we 
need to estimate the partial derivatives of / to solve equation (1). 
For this reason, it is common to integrate equation (1) over velocity 
to obtain a set of moment equations: the Jeans equations (Binney 
& Tremaine 2008). Adopting cylindrical polar coordinates ( R , <f>, z) 
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and focussing on the z-Jeans equation, we have 


Id, Id, x 1 d , 

YvM^^Tvdi^Kdz^ 


‘axial’ term: JR 



( 2 ) 


where v and <r 2 are the density and vertical velocity dispersion pro¬ 
file of a tracer population as a function of height z above the disc 
plane, cr Rz is the cross term in the velocity dispersion tensor that 
couples radial and vertical motions, and cr^ z is the cross term cou¬ 
pling vertical and axial motions. 

Integrating both sides of equation (2), we derive our key equa¬ 
tion for this paper: 

(z) = - f v( -z') K(z') - T(z') - Mz')\ dz' + (3) 

v(z) Jo Hz) 

where C is a normalisation parameter. For z = 0, we have that 
0^(0) v(0) = C, (4) 


and so C simply sets the vertical velocity dispersion at z = 0. 
(Note that a similar but less general equation was derived recently 
in Smith et al. (2012), equation (10).) We implement the constant 
C as a parameter in our model that we will ultimately marginalise 
over. The alternative would be to calculate C directly from equa¬ 
tion (3), utilising the fact that as z —» oo, v(z) —> 0, which gives us 


C = - f v(z')[^(z')-T(z')-^(z')]&'. (5) 

Jo 

However this would mean that the derived value for <t z (z) at some 
finite z would depend the properties of our mass model, tilt and 
axial terms not just in the range [0, z] but also for the range [0, oo). 

Equation (3) is numerically appealing to solve since, unlike 
many previous methods (e.g. Garbari et al. 2011, 2012), it does not 
require any numerical differentiation and is therefore rather robust 
to noise in the data. Furthermore, equation (3) is valid for any grav¬ 
ity theory and can therefore be used as a constraint on alternative 
gravity models. In this sense, we follow the early pioneering work 
of Hill (1960) who attempted to also measure K- directly without 
reference to the Poisson equation. 

To connect the vertical acceleration K z to the surface mass 
density of the disc, however, we must specify a gravitational model. 
For standard Newtonian weak field gravity, this is given by the 
Poisson equation: 


V 2 ® : 


d 2 ® 


1 dV 2 (R) 
~dd + R dR 


■■ 4nGp 


‘rotation curve’ term: H 


( 6 ) 


where V C (R) is the circular speed (rotation) curve at radius R, and 
p is the total matter density. Notice that for a flat rotation curve, 
V, = const., and the ‘rotation curve’ term R vanishes. If R does not 
vanish, then it appears as a shift (potentially as a function of height 
z) to the recovered p DM that can be corrected for at the end of the 
calculation (Garbari et al. 2012). Equation (6) can be rewritten as 

< 9 2 ® 

= 4nGp(z) eff (7) 

where 


P(z)eff = p(z) - 


1 dV 2 (R) 
4 nGR dR 


( 8 ) 


The contribution of this term to the mass density profile can be 


Table 1. Values of Oort constants. We include the F giants even though the 
errors for them are substantially larger to show that, within current uncer¬ 
tainties, the Oort constants A and B do not depend on stellar type. 


Source 

Stellar 

Type 

A 

[km s~ 

B 

1 kpc~ 1 ] 

Mignard (2000) 

K-M giants 

14.5 ± 1.0 

-11.5 ± 1.0 

Branham (2010) 

F giants 

14.85 ± 7.47 

-10.85 ±6.83 

Branham (2011) 

G giants 

14.05 ± 3.28 

-9.30 ± 2.87 


quantified via the Oort constants A and B (Binney & Tremaine 
2008): 

1 dVl(R) B 2 - A 2 

4 nGR dR 2nG ' ( 

Note that |B| < |A|, meaning that the terms in equation (9) are 
negative. Thus the effective density p(z) e ff will be higher with the 
inclusion of a non-zero rotation curve teim. If the rotation curve 
term is erroneously neglected it will be absorbed into the density 
profile, yielding an over-estimate of the DM density, assuming the 
baryonic contribution is well constrained. Taking the most precise 
values for A and B from Table 1 (Mignard 2000) yields a value 
for equation (9) of ~ 0.1 GeV/c 2 , which is roughly a third of the 
expected local DM density (e.g. Read 2014). For any accurate mea¬ 
sure of the local DM density derived from real data we will need to 
incorporate this correction, but we leave this for future work. 

Neglecting R and integrating both sides of equation (6), we 
derive the familiar result: 

£,<=> = H (10, 

where £,(z) is the surface mass density of the disc. 

The overall flow of our method is now apparent. We model 
the tracer density v(z), the mass density distibution p(z) = Pdm(z) + 
Pbaryon(z) (which gives us the K- term), the tilt term T (z), and the 
axial term tR (z). These elements each have a number of parameters, 
and so in total each model will have an A'-dimensional parameter 
space. Specifying values for each of these parameters will then give 
us a quantitative values for each of these elements, which can then 
be inserted into equation (3) to derive cr 2 (z) for that parameter space 
point. We then compare v(z) and cr 2 (z) (and cr Rz {z) as part of the tilt 
term model) to data via a^ 2 test, and then change the values of our 
parameters. 

Note that the only assumptions that goes into equation (3) are 
that the motions of stars obey the collionless Boltzmann equation, 
and that the galaxy is in dynamic equilibrium. The assumption of 
dynamic equilibrium can be negated somewhat by the presence of 
disequilibria (‘wobbling’) in the disc caused by e.g. the Sagittar¬ 
ius merger (Purcell et al. 2011; Gomez et al. 2013) or in part by 
the presence of spiral arms (Faure et al. 2014). However the im¬ 
pact of these effects on the measurement of the local DM density 
is estimated to be approximately 10% (Widrow et al. 2012; Read 
2014), less than the corrections arising from the tilt and rotation 
curve terms. 

In practice further assumptions are necessary to model the in¬ 
dividual components on the RHS of equation (3). However this 
method can accommodate almost any model for each of these ele¬ 
ments - the only strict requirement is that each element can be de¬ 
fined at the z-values corresponding to the centres of the bins used to 
calculate the tracer density and velocity dispersion. Hence we call 
this method ‘non-parametric’ - we can in principle have many more 
parameters for each model than there are data points. In the follow- 
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ing subsections we describe how we model each element for this 
particular study. In Section 2.6 we discuss model selection using 
the Bayesian evidence, which could potentially allow us to com¬ 
pare alternate models. 

As the galaxy is close to being axisymmetric in both the thin 
disc (Dehnen 1998; Hogg et al. 2005; Aumer & Binney 2009; Bovy 
et al. 2009; Pasetto et al. 2012b; Aghajani & Lindegren 2013) and 
the thick disc (Pasetto et al. 2012a; Sharma et al. 2014), the ax¬ 
ial term 3K(z) is expected to be small. For this study we assume 
complete axisymmetry, and thus take 2H(z) = 0. If the data shows 
significant non-axisymmetry in the Milky Way, the axial term could 
be modeled in a similar way to the tilt term as described below. 

Note that if the axial term were non-negligible, in the con¬ 
text of the Milky Way, this would imply the presence of significant 
spiral arms or a bar. Such features are time dependent, implying 
that if 2H(z) ± 0, we should also worry about the time dependence 
of the gravitational potential <F(r) and, therefore, also the partial 
time derivative of the distribution function df /dt (e.g. Monari et al. 
2016). It is not clear if there can exist a regime in which M(z) ± 0 
and such time dependence can be ignored. Here, we simply note 
that since empirically 3ft (z) <K 7~(z) Vz (e.g. Pasetto et al. 2012a, b), 
to leading order we can also safely neglect the time dependence of 
the distribution function, as we have already done. 


While this model is not likely realistic for the Milky Way (Flynn 
et al. 2006; McKee et al. 2015), it has been applied to observational 
data by Kuijken & Gilmore (1989a,b), and also more recently by 
Zhang et al. (2013). When applying our method to real data we 
will consider more sophisticated baryon models, but for this initial 
study where we are primarily interested in testing our methodology, 
equation 13 will suffice. 

2.4 Dark Matter parameterisation 

The simplest way to parameterise the DM density profile p DM is to 
assume it is constant with z, as done in previous work (e.g. Bahcall 
1984; Kuijken & Gilmore 1989a, b; Creze et al. 1998; Garbari et al. 
2011, 2012). This assumption works well at low z: for a spherical 
NFW halo with a scale radius of 20 kpc the midplane value is cor¬ 
rect within 10% up to a height of z ~ 3 kpc. For some analyses we 
also make this assumption, and set Pdm(z) = Pdm, const.- 

However, this assumption does not allow for the exploration 
of several interesting effects such as a flattened, oblate halo or a 
dark disc. To accommodate such phenomena we add a dark disc 
(DD) on top of the constant DM. The dark disc is described using 
the same disc model as we use for the baryons, yielding a total DM 
profile: 


2.2 Tracer density model 

To apply the Jeans equations we bin data to obtain the tracer density 
v and the velocity dispersions cr- and cr Rz . Thus at a bare minimum 
we only have to define v and cr Rz , and derive cr- via equation (3), at 
the bin centers, i.e. at a discrete set of z values. 

For this work we model the tracer density as a sum of N expo¬ 
nential discs: 


viz) 


Pdm(z) - Pdm, const. + 7Z 

4/r(_r 


KddD: 


DD 


(D 2 +Z 2 ) 1 ' 


(14) 


To ensure this disc does not simply become degenerate with the 
baryonic disc we give the scale height Ddd a higher prior range 
than that for D bn . In terms of the K- parameter of equation (3), the 
DM profile is 





2 Fz+ — 

-f . 

(11) 


I 

\ h <) 



\ 


KddZ 


k 2 + Dl D 


(15) 


where for the i' ,h disc v(0), is the tracer density at z = 0 and /z, is 
the scale height. The number of exponential disks can be increased 
until the Bayesian evidence shows no improvement (see Section 
2.6). This method avoids overfitting as each disc is smooth, while 
still giving freedom to describe more complex data. From this point 
onwards we consider only z > 0. 


2.3 Baryon parameterisation 


For this study we use a simple 2-parameter model to describe the 
baryon distribution 1 : 


Pbaryon(z) — 


1 

4 nG 


(°l + Z 2 ) 1 - 5 


( 12 ) 


The K bn parameter sets the mass of the disc, and has dimensions 
of acceleration, while D bn controls the scale height of the disc, and 
has dimensions of length. Expressed in terms of the K z parameter 
of equation (3), the baryon profile becomes 


where F = 2nGpoM,const- This disc parameterisation of the non¬ 
constant DM profile can describe both the accreted dark disc and a 
flattened DM halo, and so henceforth we refer only to a ‘dark disc’. 

Additional dark disc terms could in principle be added to 
Pdm(z) with parameters K DDji and D dd „, giving a total density pro¬ 
file of 


Pdm(z) - Pdm ,const. 


—y 

4 nG 


KdD.iiD 


2 

DDji 


(£dd,„+z 2 )‘ 


(16) 


The prior ranges on the n th dark disc parameters K DD „ and D DD n 
can be set in relation to those of the (n— l) th dark disc, e.g. requiring 
Kdd.h < K dd „_i, meaning each dark disc is less massive than the 
previous one. The number of dark discs to add could be determined 
from the data via the Bayesian evidence, with dark disc terms being 
added up until this degrades significantly. For this study however 
we limit ourselves to one dark disc term. 


2.5 Tilt term 


K- 


:,baryon : 


^bnZ 

V z2 + D b„ 


(13) 


1 In the context of the values of the prior ranges given in Table 2 below the 
gravitational constant equals G = 4.229 X 10~ 6 knr kpc M,, 1 s ' 2 . 


The tilt term from equation (2) links radial and vertical motion. The 
importance of this term has been noted previously, e.g. Kuijken & 
Gilmore (1989a); Smith et al. (2009); Biidenbender et al. (2015). 
Here we demonstrate that it is possible to deal with this term while 
remaining within our vertical, one-dimensional framework. Given 
the data quality currently available for the Milky Way we are re¬ 
quired to make a well-motivated assumption about the radial form 


MNRAS 000. 1-11 (2016) 



















A non-parametric method for measuring the local dark matter density 5 



2 [kpc] 

Figure 1. <tr z for two populations of G-type dwarf stars in the solar neigh¬ 
bourhood from the SDSS /SEGUE survey (Budenbender et al. 2015). The 
blue data points are from a younger, high metallicity population, with 
[Fe/H] = -0.07 and [a/H] = 0.11, while the data points in red are from an 
older, low metallicity population with [Fe/H] = -0.89 and [cr/H] = 0.34. 
The blue and red lines are power laws in the form of equation (20) fitted to 
the blue and red data points respectively. The older, lower metallicity stars 
probe further above the disc plane. Dashed lines indicate extrapolation from 
data. 



Figure 2. The importance of tilt, as quantified by equation (22). The blue 
and red lines correspond to tilt terms derived from the high and low metalic- 
ity population fits of Fig. 1. The deviation arising from the tilt term increases 
more rapidly with height for the low metallicity population (red line), 
which probes the high-z region most useful for determining the DM profile. 
Dashed lines indicate extrapolation from data. For this case Rq = R\ = 2.5 
kpc., and F = 267.65 km 2 kpc -2 s -2 . 


of T, but in the Gaia-era we will be able to directly measure this 
from the data. 

We first assume that the radial profiles of the tracer density 
and the (R, z)-velocity dispersion are exponentials with scale radii 
of Ro and R\, respectively: 


v(R,z) = vfe)exp(-/?//? 0 ), (17) 

<trz( R ’Z) = VR Z (.z)exp(-R/R 1 ). (18) 


With this assumption the tilt term becomes 


T(R,z) = <t r -(R,z) 


1 


1 

r! 


(19) 


If the disc were observed not to be an exponential then an alter¬ 
native model could be easily applied at this stage. Indeed, similar 
but more complex models have featured previously in the literature, 
e.g. (Binney et al. 2014). 

We then model the (R, z)-velocity dispersion as a power law at 
a given galactocentric radius R: 


<Tr z (R, z) 



( 20 ) 


We take R = R 0 , and also Ro = R i, simplifying the tilt term to a 
model decribed by the parameters A,and R 0 : 


T(R e ,z ) 


(J_\- 


1 2 

\kpcj 

Rq 

Rq Rq 


( 21 ) 


Note that we are not affected by the assumption R 0 = R t since 
these two parameters are trivially degenerate in equation 21. It 
remains the case that any observational constraints on Rq and/or 
R[ can be used as priors on equation 21, where these would con¬ 
strain the term 2/Rq. The description of cr RiZ (z) of equation (20) 
fulfils cr R '-(z = 0) = 0 by construction and fits remarkably well 
to different populations, as shown in Fig. 1. This figure presents 
(z, cr Rz ) data points for high and low metallicity populations (blue 
and red points respectively) from the SDSS /SEGUE survey as anal¬ 
ysed and presented in Budenbender et al. (2015), with a sign cor¬ 
rection applied 2 . The high metallicity sample has [Fe/H] = -0.07 
and [tr/H] = 0.11, while the low metallicity sample has [Fe/H] = 
-0.89 and [cr/H] = 0.34. Metallicity can be used as a proxy for 
age, with the high metallicity sample being younger than the older, 
low metallicity population (e.g. Loebman et al. 2011). 

Taking these data points we fit power laws models as per equa¬ 
tion (20), with parameters A = 123.99 km 2 s~ 2 , n = 1.16 for 
the high metallicity population (blue) and A = 180.08 km 2 s~ 2 , 
n = 1.44 for the low metallicity population (red). The low metallic¬ 
ity population samples further above the disc plane and populates 
the canonical thick disc of the Milky Way. Populations such as this 
are more interesting for local DM searches as they allow us to probe 
higher 2 regions where the baryon mass contribution begins to drop 
away leaving behind the DM halo. However as we go higher in the 
disc the tilt term becomes increasingly important, as illustrated by 
Figure 2, which shows the variable g(z): 


£ (z) = 


^;,DM 

r z.dm ~ 7” 


( 22 ) 


where K z DM = -2 Fz, the constant DM density term, with F = 
267.65 km 2 kpc~ 2 s -2 and Rq = Ri = 2.5 kpc, the same values as 
are used to later generate our mock data in Section 3. Compared 
to the thin disc population (blue line), the effects of the tilt term 


2 Private communication with authors. 
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Table 2. Prior ranges for parameters. Gaussian priors are described by a 
median p and a standard deviation SD. Note that v(0) and cr-(O) are the 
tracer density and velocity dispersion at z = 0, while quantities subscripted 
with 0, such as vo and SD V o, are the values derived from data in the 0 th bin, 
whose bin centre jo > 0. Tracer density v(0) has units of [stars kpc A'j,,, 
and K\)\) terms have units [km 2 kpc -1 s~ 2 ]. 


Model Parameter 

Range or Gaussian /i & SD 

Type 

Baryons 

Kbn 

p = 1500, SD = 150 

Gaussian 


An 

H = 0.18 kpc, SD = 0.02 kpc 

Gaussian 

Constant DM 

PDM 

[2,20] x 10~ 3 M 0 per 3 
[0.076,0.796] GeVcm" 3 

Linear 

Dark Disc 

A>d 

[0, 1500] 

Linear 


Odd 

[1.5, 3.5] kpc 

Linear 

Tilt term 

A 

[0, 200] km 2 s —3 

Linear 


n 

[1.0, 1.9] 

Linear 


Ro 

H — 2.5 kpc, SD - 0.5 kpc 

Gaussian 

Tracer density 

v(0) 

[0, vo + 2 x SD v o] 

Linear 


h 

[0.4, 1.4] kpc 

Linear 

Velocity disp. 

o%(0) 

[cr z , o - 2 x SDo- z , o, 
ct z q + 2 X SD^q] km s -1 

Linear 


become important at much lower z values. In short, to probe local 
DM we want to use thick disc stars probing higher-;; ranges, but the 
cost is that we must include the tilt term in our analyses. 

Ignoring the tilt term will always cause an underestimation of 
the local DM density, if all other components of the model such 
as baryons are held steady. Looking at equation (21), we note that 
R 0 < R e and A > 0, meaning the tilt term T(R a , z) is always neg¬ 
ative. Then considering equation (3) we see that to fit to cr 2 (z) in 
the absence of the tilt term T (z'), the K,{z') term, a negative term, 
is forced to become less negative in order to compensate. This re¬ 
quires a lower mass density, which if the baryon density profile is 
held constant, manifests itself as a decrease in the DM mass den¬ 
sity. 


2.6 Statistical analysis and MultlNest 


The model outlined above gives us an (V-dimensional parameter 
space. To explore this parameter space and derive limits on the var¬ 
ious observables we adopt nested sampling (Skilling 2006) as im¬ 
plemented in the publicly available MultiNest code (Feroz & Hob¬ 
son 2008; Feroz et al. 2009; Feroz et al. 2013). MultiNest is a tool 
for Bayesian inference and parameter estimation. Bayes Theorem 
is 


P(6\D, M) 


P(D\9,M)P(9\M) 

P(D\M) 


(23) 


where At is the given model, 9 is the set of parameters for that 
model, and D is the data. The left hand side, P(6\D, M) is known as 
the posterior , while the three terms on the right are the likelihood 
P(D\9, M) = -L(6), the prior P(9\M), and the Bayesian evidence 
P(D\M). 

The Bayesian evidence, a.k.a. the marginal or model likeli¬ 
hood, is a measure of how well the model performs given the data, 
and can be expressed as 

Z = P(D\M)= I P(D\9, M)P(9\M)d9. (24) 


The performance of two different models given the same data can 
be compared using the Bayes factor. 


P(D\Mo ) 
P{D\Mi) 


(25) 


Assuming Gaussian errors it is possible to derive an empirical 
scale relating the Bayes factor to the strength of evidence for one 
model over another, as done in (Trotta 2008). There, a | In B 0 il value 
of less than 1 is considered inconclusive, while values of 1.0, 2.5 
and 5.0 are considered to give weak, moderate and strong evidence, 
respectively. 

MultiNest takes a given prior probability distribution and 
likelihood function and calculates the posterior distribution and 
Bayesian evidence. Our likelihood function is based on the^ 2 dis¬ 
tribution: 

£(0) = ex pf-f)- (26) 


The value of x 2 is simply 


2 2 2 2 

T -Tv + X a 2 +To- ff -' 

where 


w 

II 

(Tdata,/ Vmodelj)^ 

SD 2 

J 


= Z 

(^"z,dataj CT 7 , model J ) 

SD 2 , 

j 


A = 2 

j 

(^D?z,dataJ — model, j) 

SD l Rl ., 


(27) 


(28) 

(29) 

(30) 


Note that if the reconstruction model does not contain a tilt term 
(e.g. T = 0) then ;y 2 is set to zero. 

Table 2 shows the prior ranges used for our analyses. We de¬ 
rive credible regions (CRs) for the DM density parameters by tak¬ 
ing its posterior distribution and marginalising over the remain¬ 
ing parameters. As outlined above our model has several compo¬ 
nents that can be turned on or off, such as the dark disc. Using the 
Bayesian evidence as calculated by MultiNest it is potentially pos¬ 
sible to perform model comparison to determine which reconstruc¬ 
tion model best fits the data. This idea will be explored in greater 
depth in subsequent studies. 


3 MOCK DATA SETS 

In this paper we apply our method only to mock data in order to 
hone and verify it. This mock data is ‘as good as it gets’, in the sense 
that it has no measurement errors, nor observational biases added 
to it, and is drawn from relatively simple, known distribution func¬ 
tions. While dynamically realistic ‘N-body’ mocks are preferred, 
it has already been shown that ID Jeans analyses are robust to the 
presence of local non-axisymmetric structure in the disc (Garbari 
et al. 2011 ). Furthermore, N-body mocks are expensive; even state- 
of-the-art simulations do not approach the local sampling expected 
from Gaia. Finally, the expense of building well sampled N-body 
mocks makes it challenging to explore a large parameter space of 
models, including models with and without tilt, or with and without 
a dark disc. For these reasons, we focus here on simpler mock data, 
similarly to the Read (2014) review. We will consider more dynam¬ 
ically realistic mocks, and mocks that give a faithful representation 
of the expected Gaia data uncertainties in future work. 

We generate here a variety of mock data sets as outlined in 
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Table 3, with different samplings (10 4 , 10 s , and 10 6 tracer stars), 
with and without a tilt term, and also with either no dark disc, a 
dark disc (dd), or a more massive ‘big dark disc’ (bdd). We assume 
the axial term JK and the rotation curve term 7? are zero. Our simple 
baryonic disc model is set up to mimic the real Milky Way, with 
a scale height parameter of D bll = 0.18kpc and a surface density 
of Z bn = 55.53 M G pc' 2 , similar to those measured near the solar 
position (Flynn et al. 2006; Read 2014; McKee et al. 2015). The 
value of the F parameter (see equation 15) corresponds to a DM 
density of Pdm, const = lOx 10' 3 M 0 pc' 3 = 0.38 GeV cm' 3 . For each 
scenario we generate 20 mock data sets, allowing us to explore the 
effects of poisson noise over a range of realisations. 

The mock data consist of a list of stars each with three pieces 
of data: the position z, the vertical velocity v z , and the product of the 
vertical and radial velocities, v R v z . The first element, z, is generated 
by drawing randomly from an exponential tracer distribution with 
scale height h: 

v(z) = exp (31) 

This gives us our list of stellar positions. 

To derive the velocites for the mock catalogue we first must 
define a mass model and tilt. Taking the same parameterisations as 
described in Sections 2.3, 2.4, and 2.5 for baryons, DM, and tilt, 
respectively, we set their parameters as per Table 3. This allows 
us to calculate K- and T(z). Using equation (5) and its associated 
assumptions we can calculate C. We then use equation (3) to derive 
cr,(z). For each star we take its position z', find the value of o\(z'), 
and then draw a velocity from a Gaussian centred on v z = 0 and 
with variance cr 2 (z'). 

To generate v R v z mock data, when necessary, we take the A, 
n , and R 0 parameters listed in Table 3 and calculate the cr R -(z) pro¬ 
file via equation (20). Taking each star’s position z', we calculate 
cr Rz (z'), and draw a value of v R v z from a Gaussian centered on 
v R v z = 0 and with variance of cr 2 R Xz')- 


4 RESULTS 

Here we present the results of our scans. We first investigate how 
the precision of the reconstruction varies with different numbers of 
tracer stars. We then look at the effects of the tilt term and the dark 
disc. For certain parts of these analyses the model used to generate 
the mock data and the model used to reconstruct the mock data are 
not the same - this is done to investigate the effects of ignoring 
terms (as in the case of the tilt term) or using incorrect models. 

In the following figures we plot the marginalised posterior for 
the DM density PdmU)> showing 68%, 95%, and 99.7% CRs in 
dark, medium and light shading, respectively. The mock data pro¬ 
file is shown by the solid black line, while the median of the pos¬ 
terior distribution is shown by a solid line in the same colour as 
the CR. Binning of stars to calculate v(z), cr-(z), and cr Rz {z) is per¬ 
formed such that each bin contains an equal number of stars. For 
this study we use 20 bins, and in Appendix A we briefly explore the 
effects of changing the number of bins used. For plots with constant 
DM density in both mock data set and reconstruction model we plot 
all 20 mock data sets together, while for non-constant DM density 
in either mock or reconstruction we show one representative figure 
in this section, and show the full set of figures in Appendix B (Figs. 
B1 to B6). 

Our method and code is set up to describe the tracer density 
as a sum of exponentials. To determine the necessary number of 
exponentials required to describe the data we can use the Bayes 


Table 3. Mock data parameters. _X is the number of stars sampled, e.g. 
10 4 , 10 5 , 10 6 , and _M is the mock number, ranging from 0 to 19. Empty 
spaces indicate that a certain mock data set does not include that ele¬ 
ment. The baryon model corresponds to a baryonic surface density of 
Z bn = 55.53 M q pc' 2 , while the F parameter corresponds to a constant DM 
density of PDM.const = 10 X 10' 3 M e pc' 3 = 0.38 GeV cm' 3 . K bn and K^o 
terms have units [km 2 kpc' 1 s' 2 ], while F has units [km 2 kpc' 2 s' 2 ]. 



Tracer 

density 

1 n n 




Potential 



^bn 

1 jUU 


Dbn LkpcJ 

F - 

- 267.65 - 


Dark Disc 

Kom 
Z>dd [kpc] 

300 

2.5 

900 

2.5 


900 

2.5 

Tilt Term 

A [km 2 s' 2 ] 



180.08 

180.08 


n 



1.44 

1.44 


Ro [kpc] 



2.5 

2.5 


factor as described earlier in Section 2.6. Test reconstructions give 
a Bayes factor of 1.5 when comparing one exponential to two, and 
3.3 when comparing one exponential to three, with one exponential 
favoured in both cases. This is to be expected given that the mock 
data is generated using a single exponential. Given this result we 
use a single exponential for all subsequent reconstructions. 


4.1 Sampling 

Fig. 3 shows reconstructions of the local DM density for using 
varying numbers of tracer stars. The mock data sets used here are 
the most simple case. thick_X_M as described in Table 3, and has 
no dark disc or tilt added. As expected the credible regions shrink as 
the number of tracer stars is increased from 10 4 stars up to 10 s stars. 
The SDSS survey has sampling on the order of 10 4 stars (Zhang 
et al. 2013; Biidenbender et al. 2015) while data from Gaia will 
give upwards of 10 6 tracer stars. 


4.2 Tilt 

In Fig. 4 we explore the effects of the tilt term. The left hand set 
of CRs in Fig. 4 are the same as the right hand set of CRs from 
Fig. 3, with no tilt term in the mock data or reconstruction. In the 
centre set of Fig. 4 however, the mock data thick_tilt_le6_M 
has the tilt term turned on, but the analyses are performed with the 
tilt term set to zero. This illustrates the danger of ignoring tilt as 
discussed earlier in Section 2.5: the reconstructions return narrow 
credible regions, but as expected they systematically underestimate 
the value of p DM , with the true p DM lying outside even the 99.7% 
CRs. This underestimation however is remedied by turning on the 
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I posterior median H Pdm, mock inside 68% CR _ 
68% CR I Pdm, mock inside 95% CR 

95% CR BpDM.mock inside 99.7% CR - 
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1 
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a 


o 
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Figure 3. Determination of the DM density profile for varying numbers of tracer stars. These plots show marginalised posteriors for pdm(z) = PDM,const for the 
20 mock data sets generated for each sampling level. Dark, medium, and light shading show the 68%, 95%, and 99.7% credible regions (CRs), respectively. 
Green, purple, and orange colouring indicates that the 68%, 95%, or 99.7% CR respectively contains the correct answer. The median value of each posterior 
is shown by a solid line in green, purple, or orange, while the DM density value used to generate the mock data is shown as a solid black line across the entire 
plot. The mock data and reconstruction models used contain no tilt term or dark disc terms. As the number of tracer stars is increased from 10 4 to 10 6 the 
credible regions for Pdm, const naturally shrink around the mock data value. 


tilt term in the model, as shown in the right hand set CRs in Fig. 4, 
where the correct DM density always at least within the 95% cred¬ 
ible region. The inclusion of extra parameters to describe the tilt 
term does however increase the size of the CRs. 

Just as our determination of pdm(z) is dependent on the tilt 
term, the tilt term is in turn dependent on its input parameters, A, 
n, and R 0 = Ri. While we have been able to fit for A and n using 
G-type dwarf data (section 2.5), for Rq we have taken the canonical 
value of R 0 = 2.5 ± 0.5 kpc from (Binney & Tremaine 2008), and 
further made the assumption that ct r -(R , z) has the same radial pro¬ 
file as v(R,z ) (i.e. Rq = Ri). When using only a single population, 
determination of R 0 and R [ will be important, as illustrated in Fig¬ 
ure 5. This figure shows the 2D marginalized posterior for pnM.const 
and Rq for a reconstruction of mock data set thick_tilt_le6_0 
with a model including tilt, and demonstrates the degeneracy that 
exists between pDM.const an d Ro = Ri. If using multiple tracer popu¬ 
lation, each will have different Rg and/or Rj , but will all have their 
motions dictated by the same potential. This will help us break the 
degeneracy between R 0 , R lt and p DM . Further, with the advent of 
Gaia data we will be able to directly measure the radial profile of 
cr Rz and v(R, z) for a given set of tracer stars. 

4.3 Dark Disc 

Fig. 6 shows mock data set thick_le6_0-19, thick_dd_le6_0, 
and thick_bdd_le6_@ reconstructed using models with and with¬ 
out a dark disc component. The top left panel shows the same CRs 
as seen in the left hand set o Fig. 4 and the right hand set of Fig. 3. 

The left column of Fig. 6 shows the reconstruction of mock 
data sets with a constant DM density profile; mocks 0-19 in top 
left and mock 0 in bottom left. The reconstruction in the top panel 
uses a model with a constant DM density, while the bottom panel 
uses a model with an additional dark disc component. The dark 
disc reconstruction exhibits a disc structure even though the correct 
answer is constant p DM . This is likely due to a bias in the hyper¬ 
volume set by the priors - the prior range on the dark disc parame¬ 
ters goes between no dark disc (K DD = 0) and a maximal dark disc 
(K DD = 1500), and thus the bulk of the parameter space features at 


least some dark disc. There is no ‘negative dark disc’ to counteract 
this effect and push the mean of the prior range back to no dark 
disc. 

In the centre and right columns of Fig. 6 we reconstruct mock 
data sets with a dark discs of different masses: thick_dd_le6_0 
and thick_bdd_le6_0 (the ‘big dark disc’). A constant DM den¬ 
sity reconstruction (top row) is able to contain the thick_dd_le6 
dark disc within the 95% credible region almost to the last bin, but 
fails to contain the big dark disc beyond z = 1.3 kpc. Adding a dark 
disc term allows the reconstruction to fit to the mock data DM pro¬ 
file nicely, as shown in the bottom centre and bottom right panels 
of Fig. 6. 

When working with real data however, we will not have the 
luxury of knowing the DM profile underlying the data. The black 
mock data model line will not be there. In subsequent studies we 
will explore the potential use of the Bayesian evidence, as calcu¬ 
lated by MultiNest , to determine if a dark disc is justified by ob¬ 
servational data. 


4.4 Tilt and Dark Disc 

Here we combine the two elements discussed in previous sections, 
the tilt term and the dark disc. Fig. 7 shows reconstructions of the 
mock data set thick_bdd_tilt_le6_0. In the top panel the re¬ 
construction model contains neither dark disc nor tilt term. Again 
we see the same effects as we did previously. The missing tilt term 
yields a consistent underestimation of the DM density, and the con¬ 
stant DM density envelope is too narrow to encompass the density 
range of the dark disc. The consistent underestimation can reme¬ 
died by adding a tilt term to the reconstruction model, as shown 
in the middle panel, where the reconstruction model has a tilt term 
included. Both problems can be resolved in tandem by using a re¬ 
construction model with both tilt and dark disc terms, as shown in 
the bottom panel of Fig. 7. 
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Figure 4. Determination of the DM density profile using 10 6 tracer stars and exploring the effects of including or neglecting the tilt term in the mock data 
sets and reconstruction models. These plots show marginalised posteriors for pdm(z) = Pdm, const for the 20 mock data sets generated for each tilt scenario. 
Dark, medium, and light shading show the 68%, 95%, and 99.7% credible regions (CRs) respectively. Green, purple, and orange colouring indicates that the 
68%, 95%, or 99.7% CR respectively contains the correct answer, while pink colouring indicated that the correct answer lies outside even the 99.7% CR. The 
median value of each posterior is shown by a solid line in green, purple, or orange, while the DM density value used to generate the mock data is shown as a 
solid black line across the entire plot. The left set of mocks shows the same result as the right set of Fig. 3, with no tilt term in mock data or reconstruction 
model. The centre set has a tilt term in the mock data but not in the reconstruction model, yielding a systemic underestimation of pdm- The right set has a tilt 
term in both mock data and reconstruction model, demonstrating that our method can successfully deal with the tilt term. 
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Figure 5. 2D marginalized posterior of the constant DM density pdm(z) = 
PDM,const and the Rq parameter of the tilt term (equation 19), generated from 
mock data set thick_tilt_le6_8 reconstructed using a model with tilt. 
Contours show the 68%, 95%, and 99.7% CRs. Marginalization and plotting 
performed using Barrett (Sebastian Liem, private communication). 


5 CONCLUSIONS 

In this paper we have presented a new method of determining the 
vertical DM density profile, and thus the local DM density. The key 
equation of this method, equation (3), depends only on the assump¬ 
tion of dynamical equilibrium. In practice further assumptions are 
made to describe the components going in to equation (3), however 
the scope of possible models is wide, and in principle model selec¬ 


tion using the Bayesian evidence can be used to determine the best 
model. We leave exploration of this last aspect to future work. 

The importance of the tilt term has been previously noted Kui- 
jken & Gilmore (1989a); Smith et al. (2009); Budenbender et al. 
(2015), and to derive an accurate value for the local DM density it 
is vital that our method be able to deal with this term. The bary- 
onic contribution to the galactic density profile drops rapidly as z 
increases, leaving DM as the dominant component. Thus to probe 
the DM density profile sampling thick disc stars that travel higher 
above the plane is preferable. For these populations the tilt term is 
more important, as illustrated in Fig. 2. Failure to include the tilt 
term in the analysis leads to a systematic underestimation of the 
local DM density, as explained in Section 2.5 and demonstrated in 
Section 4.2. 

One of the novel aspects of our method is that it can deal 
with the tilt term while remaining within the confines of the one¬ 
dimensional ^-direction Jeans equation, which can be seen in Fig. 4. 
With only the data currently available for the Milky Way, this re¬ 
quires several well motivated assumptions, as described in Section 
2.5. However, with data from Gaia we will be able to directly mea¬ 
sure the radial profile of the tilt and tracer density, removing the 
need for such assumptions. While for this paper we have disre¬ 
garded the rotation curve term (equation 6), we note that an accu¬ 
rate determination of this will be necessary for the implementation 
of this or other ^-direction methods to real data. 

Non-spherical DM density profiles, such as oblate halos or 
accreted dark discs, can also be fitted using our method by incor¬ 
porating a dark disc term into the reconstruction model, which is 
shown in Fig. 6. Our method can also reconstruct mock data sets 
containing both a tilt term and a dark disc, as shown in Fig. 7. 
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Figure 6. Determination of the DM density profile using 10 6 tracer stars and exploring the effects of including or neglecting a dark disc in the mock data sets 
and reconstruction models. For comparison the top left panel shows the same reconstructions as the right hand set from Fig. 3, i.e. basic mock data sets with 
no DD, reconstructed with a constant DM density. The remaining panels each show one representative mock and reconstruction, with the full set of mocks 
and reconstructions given in Appendix B. These panels show marginalised posteriors forpouRz). with dark, medium, and light shading indicating the 68%, 
95%, and 99.7% CRs respectively. The median value of the posterior is shown by the solid blue line, while the DM density profile used to generate the mock 
data is shown as a solid black line. From left to right the columns show reconstructions of mock data sets containing no dark disc (pdm, const only), a dark disc 
(DD), or a ‘big’ dark disc (BDD). The top centre and top right panels show the determination of a constant DM density from the DD and BDD mock data sets 
respectively. The bottom row shows the reconstruction of each of the mock data sets using a model with a constant DM term and a dark disc term. Light grey 
vertical lines indicate the bin centres. 
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Figure 7. Determination of the DM density profile using 10 s tracer stars 
with a combination of a tilt term and a dark disc in the mock data 
(thick_bdd_tilt_le6_8) and using a variety of reconstruction models. 
These plots show marginalised posteriors for pdm(z), with dark, medium, 
and light shading indicating the 68%, 95%, and 99.7% CRs respectively. 
The median value of the posterior is shown as the solid blue line, while the 
DM density profile used to generate the mock data is shown as a solid black 
line. 
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APPENDIX A: VARIATION DUE TO THE NUMBER OF 
BINS 

Figure Al illustrates the effects of changing the number of bins 
used for this analysis. There we plot the 68%, 95%, or 99.7% CRs 
for 20 mock data sets (thick_lE6_0-19, no tilt and no DD), and 
vary the number of bins from five to 30, in increments of five. For 
only five bins (top left set), the true answer for Pdm, const is outside 
the 99.7% CR for all but two of the mocks, for which the true an¬ 
swer is within only the 99.7% CR. The systematic underestimation 
is due to an over estimation of the baryonic disk. The fifth bin in 
this scheme covers a range from 1.2 kpc to 2.4 kpc and has its bin 
centre at z = 1.6kpc, so it is unsurprising that such a low number 
of bins fails to correctly reconstruct the DM profile, which is a sub¬ 
dominant component only becomes apparent at higher z. Increasing 
the number of bins to 10, 15, and then to 20 improves the results. 
The gains from increasing from 20 to 25 and 30 bins is very slight. 


APPENDIX B: DARK DISC RECONSTRUCTIONS WITH 
MULTIPLE MOCKS 

Here we show the results of reconstructing all mocks data sets gen¬ 
erated using a constant DM density only (thick_X_M), constant 
DM and a regular dark disc (thick_dd_X_M), or constant DM and 
a ‘big’ dark disc (thick_bdd_X_M); see Table 3 for details. The 
reconstruction models either have a constant DM density only, or a 
constant DM density plus a dark disc. 

This paper has been typeset from a TpX/IATj:X file prepared by the author. 


MNRAS 000. 1-11 (2016) 


A non-parametric method for measuring the local dark matter density 13 


"I 

I 


s 

§ 


I 


posterior median I PDM.mock inside 68% CR 
68% CR HpDM.mock inside 95% CR 

95% CR Bp DM ,ra OC k inside 99.7% CR 

99.7% CR I PDM ,mock outside 99.7% CR 




0.4 ^ 

O 




5 bins 


10 bins 


15 bins 


a 

L 


2 

§ 


|l l ll l ll|llH l t||l|l |I I II I I | |I|H I| I| 1 |I |l l H l |l|lll |l *l|l|l 


> 
: CD 

o 


20 bins 


25 bins 


30 bins 


Figure Al. Exploring the effects of binning on the determination of pdm- These plots show marginalised posteriors for Pdm(z) = Pdm, const f° r the 20 mock 
data sets thick_lE6_Q-19, reconstructed using 5, 10, 15, 20, 25, and 30 bins. Dark, medium, and light shading show the 68%, 95%, and 99.7% credible 
regions (CRs) respectively. Green, purple, and orange colouring indicates that the 68%, 95%, or 99.7% CR respectively contains the correct answer, while 
pink colouring indicated that the correct answer lies outside even the 99.7% CR.. The median value of each posterior is shown by a solid line in green, purple, 
or orange, while the DM density value used to generate the mock data is shown as a solid black line across the entire plot. 
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Figure Bl. Marginalized posteriors of pdm(^) for mock data sets thick_lE6_Q - 19 with pdm, const only (no DD component) reconstructed using a model with 
Pdm, const plus a DD. Dark, medium, and light shading indicate the 68%, 95%, and 99.7% CRs respectively. The median value of the posterior is shown as the 
solid blue line, while the DM density profile used to generate the mock data is shown as a solid black line. 
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Figure B2. Marginalized posteriors of pdm(z) for mock data sets thick_dd_lE6_©-19 with Pdm, const plus a DD component, reconstructed using a model 
with pdm, const only. Dark, medium, and light shading indicate the 68%, 95%, and 99.7% CRs respectively. The median value of the posterior is shown as the 
solid blue line, while the DM density profile used to generate the mock data is shown as a solid black line. 
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Figure B3. Marginalized posteriors of pdm(z) for mock data sets thick_dd_lE6_©-19 with Pdm, const plus a DD component, reconstructed using a model 
with Pdm, const and a DD component. Dark, medium, and light shading indicate the 68%, 95%, and 99.7% CRs respectively. The median value of the posterior 
is shown as the solid blue line, while the DM density profile used to generate the mock data is shown as a solid black line. 
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Figure B4. Marginalized posteriors of pdm(z) for mock data sets thick_bdd_lE6_0-19 with Pdm, const plus a DD component, reconstructed using a model 
with pdm, const only. Dark, medium, and light shading indicate the 68%, 95%, and 99.7% CRs respectively. The median value of the posterior is shown as the 
solid blue line, while the DM density profile used to generate the mock data is shown as a solid black line. 
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Figure B5. Marginalized posteriors of pdm(z) for mock data sets thick_bdd_lE6_0-19 with pdm, const plus a ‘big’ DD component, reconstructed using a 
model with pDM,const and a DD component. Dark, medium, and light shading indicate the 68%, 95%, and 99.7% CRs respectively. The median value of the 
posterior is shown as the solid blue line, while the DM density profile used to generate the mock data is shown as a solid black line. 
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Figure B6. Marginalized posteriors of pdm(z) for mock data sets thick_bdd_tilt_X_Q-19 with Pdm, const plus a ‘big’ DD component and including the 
elfects of tilt, reconstructed using a model with pdm, const, a DD component, and incorporating tilt. Dark, medium, and light shading indicate the 68%, 95%, 
and 99.7% CRs respectively. The median value of the posterior is shown as the solid blue line, while the DM density profile used to generate the mock data is 
shown as a solid black line. 
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